﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace new_problem
{
    class MiscUtils
    {
        public static double Norm(double[,] A)
        {
            double sum = 0.0;
            for (int i = 0; i < DataGen.N + 1; i++)
                for (int j = 0; j < DataGen.N + 1; j++)
                    sum += A[i, j] * A[i, j];
            return Math.Sqrt(sum);
        }

        public static Double[,] SubstractMatrix(Double[,] A, Double[,] B)
        {
            Double[,] sub = new double[DataGen.N + 1, DataGen.N + 1];
            for (int i = 0; i < DataGen.N + 1; i++)
                for (int j = 0; j < DataGen.N + 1; j++)
                    sub[i, j] = A[i, j] - B[i, j];
            return sub;
        }

        public static Double[,] MultiplyMatrix(Double[,] A, Double[,] B)
        {
            Double[,] mult = new double[DataGen.N + 1, DataGen.N + 1];
            for (int i = 0; i < DataGen.N + 1; i++)
                for (int j = 0; j < DataGen.N + 1; j++)
                    for (int k = 0; k < DataGen.N + 1; k++)
                        mult[i, j] += A[i, k] * B[k, j];
            return mult;
        }

        public static Double[] MultiplyMatrixVector(Double[,] A, Double[] x)
        {
            Double[] mult = new double[DataGen.N + 1];
            for (int i = 0; i < DataGen.N + 1; i++)
            {
                double temp = 0.0;
                for (int j = 0; j < DataGen.N + 1; j++)
                    temp += A[i, j] * x[j];
                mult[i] = temp;
            }
            return mult;
        }

        // регуляризация по Тихонову - см. метод Левенберга-Марквардта
        public static Double[,] TikhonovRegularization(Double[,] A, double alpha)
        {
            for (int i = 0; i < DataGen.N + 1; i++)
                A[i, i] += alpha;
            return A;
        }

        // для (мин. ошибки) 
        public static double CalcNumerator(Double[,] A, Double[,] g)
        {
            return Math.Pow(Norm(SubstractMatrix(A, g)), 2);
            //  ||A-g||^2
        }

        // для (наискор. спуска) 
        public static double CalcDenominator(Double[,] dA, double[,] S)
        {
            return Math.Pow(Norm(MultiplyMatrix(dA, S)), 2);
            // ||dA*S||^2
        }

        // для градиентных методов (мин. ошибки; наискор. спуска)
        public static Double[,] CalcS(Double[,] A, Double[,] dA, Double[,] g)
        {
            // dA* (A(u)-g)
            return MultiplyMatrix(dA, SubstractMatrix(g, A));
        }

        public static void InvertMatrix(Double[,] A)
        {
            int[] row = new int[DataGen.N];
            int[] col = new int[DataGen.N];
            Double[] temp = new Double[DataGen.N];
            int hold, I_pivot, J_pivot;
            double pivot, abs_pivot;

            // установиим row и column как вектор изменений.
            for (int k = 0; k < DataGen.N ; k++)
            {
                row[k] = k;
                col[k] = k;
            }
            // начало главного цикла
            for (int k = 0; k < DataGen.N ; k++)
            {
                // найдем наибольший элемент для основы
                pivot = A[row[k], col[k]];
                I_pivot = k;
                J_pivot = k;
                for (int i = k; i < DataGen.N ; i++)
                {
                    for (int j = k; j < DataGen.N ; j++)
                    {
                        abs_pivot = Math.Abs(pivot);
                        if (Math.Abs(A[row[i], col[j]]) > abs_pivot)
                        {
                            I_pivot = i;
                            J_pivot = j;
                            pivot = A[row[i], col[j]];
                        }
                    }
                }
                if (Math.Abs(pivot) < 1.0E-10)
                {
                    throw new Exception("Matrix is singular !");
                    return;
                }
                //перестановка к-ой строки и к-ого столбца с стобцом и строкой, содержащий основной элемент(pivot основу)
                hold = row[k];
                row[k] = row[I_pivot];
                row[I_pivot] = hold;
                hold = col[k];
                col[k] = col[J_pivot];
                col[J_pivot] = hold;
                // k-ую строку с учетом перестановок делим на основной элемент
                A[row[k], col[k]] = 1.0 / pivot;
                for (int j = 0; j < DataGen.N; j++)
                {
                    if (j != k)
                    {
                        A[row[k], col[j]] = A[row[k], col[j]] * A[row[k], col[k]];
                    }
                }
                // внутренний цикл
                for (int i = 0; i < DataGen.N; i++)
                {
                    if (k != i)
                    {
                        for (int j = 0; j < DataGen.N; j++)
                        {
                            if (k != j)
                            {
                                A[row[i], col[j]] = A[row[i], col[j]] - A[row[i], col[k]] * A[row[k], col[j]];
                            }
                        }
                        A[row[i], col[k]] = -A[row[i], col[k]] * A[row[k], col[k]];
                    }
                }
            }
            // конец главного цикла

            // переставляем назад rows
            for (int j = 0; j < DataGen.N; j++)
            {
                for (int i = 0; i < DataGen.N; i++)
                {
                    temp[col[i]] = A[row[i], j];
                }
                for (int i = 0; i < DataGen.N; i++)
                {
                    A[i, j] = temp[i];
                }
            }
            // переставляем назад columns
            for (int i = 0; i < DataGen.N; i++)
            {
                for (int j = 0; j < DataGen.N; j++)
                {
                    temp[row[j]] = A[i, col[j]];
                }
                for (int j = 0; j < DataGen.N; j++)
                {
                    A[i, j] = temp[j];
                }
            }
        }

        public static Double[,] TransposeMatrix(Double[,] A)
        {
            double tmp;
            for (int i = 0; i < DataGen.N + 1; i++)
            {
                for (int j = 0; j < i; j++)
                {
                    tmp = A[i, j];
                    A[i, j] = A[j, i];
                    A[j, i] = tmp;
                }
            }
            return A;
        }

    }
}
